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I. INTRODUCTION 



A keystone procedure to obtain macroscopic thermodynamics quantities (e.g., energy, 
specific heat, magnetization, phase transition points, etc) of statistical systems is to per- 
form appropriate averages over their microscopic configurations. In practice, however, such 
systems usually have a prohibitive number of states for a full covering. Therefore, ap- 
proaches relying on proper representative samplings must be considered and so Monte Carlo 
tools become fundamental for calculations. By a proper sampling we mean that for a given 
instance a method should satisfactorily: (i) represent the way the system actually evolves 
throughout the different microstates (among the whole set S of microstates in the system); 
and (ii) generate a set Q of visited microstates that indeed gives a good picture of all the 
relevant microstates which describe the problem at that particular situation. 

Within this framework, an important issue is to know under what conditions the above 
criteria are fulfilled. For example, biased values for physical quantities may arise when the 
system displays local free-energy minima and the dynamics used to generate the microscopic 
configurations either is not able to cross such barriers or it does so, but only after too 
long times. Consequently, we have broken ergodicity for finite (even large) simulations^, 
leading to metastability and thus to poor estimates for the system properties due to a non- 
representative Q. Metastability and broken ergodicity appear in several problems like; spin- 
glasses; protein folding, biomolecules; and random search, to name just a few^. Moreover, 
they are not restricted only to complex systems, also being present in simpler contexts 
like in lattice-gas models displaying first-order phase transitions 4 - - -. As noted, in such case 
the sampling dynamics may present difficulties to cross the energetic barriers. Then, the 
system can develop hysteresis by passing back and forth the phase frontiers as we change 
the parameter control 4 -. 

Different alternative ideas have been considered to overcome^ or even circumvent^ en- 
tropic barriers, thus restoring the ergodic behavior. In particular, enhanced sampling algo- 
rithms, such as parallel tempering (PT)£ - — - also known as multiple replica exchange - and 
simulated tempering (ST)^" 1 ^, have recently attracted a lot of attention, specially due to 
their simplicity and generality compared to other Monte Carlo algorithms 4,5 . Briefly, in the 
PT method, microscopic configurations in higher temperatures are used to assure an ergodic 
free walk in lower temperatures: one simulates replicas of the same system at distinct T's, 
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allowing the exchange of temperature between the replicas. For the ST, on the other hand, 
an unique replica is considered, however, the system occasionally undergoes temperature 
changes along its evolution. 

Given the different tempering implementation in the two approaches, a natural question 
is how they compare to each other— ~— . For example, the rate of temperatures switching is 
higher for the ST— ~— . So, usually one could expect a larger number of distinct phase space 
regions visited when using the ST, thus a possible advantage over the PT. But as we discuss 
in Section II. C, near phase transition conditions this is not always the case. Therefore, it 
still an open query if indeed one method is systematically superior in all situations. With 
the above in mind, here we compare the PT and ST efficiencies when applied to phase 
transitions, specially to the first order case. 

A short comment regarding the comparison between the PT and ST for first order phase 
transitions is in order. In principle, for a true first order transition, i.e., for systems in the 
thermodynamic limit, the energy descontinuous gap would lead to a small probability of 
accepting exchanges between the PT replicas^. But in concrete calculations, one is always 
dealing with finite sizes L, where the actual thermodynamics properties are described by 
continuous functions. Also, these functions are smooth and tend to the correct asymptotic 
behavior (for L — > oo) only if the state space is properly sampled^ 7 , what has been shown 
to be the case for the PT— . Thus, in practice the above mentioned difficulty for the PT is 
not an issue and the method is indeed an appropriate tool to study first order transitions, 
as discussed and exemplified in different works^&i 7 -. Hence, the PT and ST (this latter 
rarely considered in such regime, few exceptions being Refs.— ) can be analyzed at the same 
footing. So, possible convergence differences can be associated just to the way the algorithms 
generate the sets Q, and not to the approaches eventual instrinsic distinctions (recall that 
conceptually they are similar—). 

In this contribution we first revisit the simplest Ising spin model displaying a well un- 
derstood second order phase transition. This is an instructive example because in a recent 
work—, it has been shown that through an improved version of the ST, the frequency of 
successful exchanges (measured in terms of transition decay rates) is higher for the ST than 
for the PT method. However, the comparison was not carried near the critical temperature. 
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By analyzing time correlation functions, defined as 

C w (r) = (w(t) - w) (w(t + t)-w), (1) 

for w relevant thermodynamic quantities (like energy and magnetization) of mean w and ( ) 
denoting time averages, one no longer gets a better performance of the ST around T c . In 
fact, we find that the PT leads to faster decaying C's. 

Then, we move to the main focus of this contribution: the harder situation of strong 
first-order phase transitions, where the use of one-flip algorithms like Metropolis often gives 
rise to poor numerical simulations. As the specific case study, we consider the lattice gas 
model with vacancies (a spin-1 model in the magnetic systems jargon) 20 . This class of 
problems has been extensively studied under different alternative methods^-! 2 ^ 2 .. Hence, 
the many available results can help to benchmark those obtained from the PT and ST. 
We show that although both, PT and ST, lead to equivalent good results in the limit of 
long simulations, the PT displays a faster convergence towards stationarity. Moreover, for 
the PT, the tunneling between different phases at the coexistence is more frequent and the 
generated microscopic configurations uncorrelate faster. 

The work is organized as the following. In Sec. II we review the PT and ST methods, 
discussing distinct implementations. We also give reasons why the PT may outperform ST 
near phase transition conditions. In Sec. Ill we consider a spin system displaying a second- 
order phase transition. The lattice-gas model and its comparative study with the PT and 
ST methods - addressing a first order phase transition - are presented in Sec. IV. Finally, 
in Sec. V we draw our last remarks and the conclusion. 

II. THE PT AND ST SAMPLING ALGORITHMS 

The central idea behind a tempering enhanced sampling algorithm is try to guarantee 
ergodicity by means of appropriate temperature changes during the simulations, thus allow- 
ing efficient and uniform visits to a fragmented multiple regions phase space 19 . Suppose we 
shall study a system at a given To. We assume T\ = To and define a set of N distinct tem- 
peratures Ti < T 2 < . . . < T/v, with AT = Tjy — T\. There are different ways to implement 
tempering 2 ^, two important ones being the PT and ST, which we describe next. 
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A. Parallel Tempering 



The PT approach combines a standard algorithm (e.g., Metropolis) with the simultaneous 
evolution of iV copies of the system (each at a different T n ), occasionally allowing the replicas 
to exchange their temperatures. Fixing relevant parameters, the method is implemented by 
first running M eq times (to assure equilibration of all the N copies) a two parts procedure, 
(a) and (b), discussed below. After that, for each (a)-(b) composite MC step (repeated M a b 
times) we calculate the thermodynamics quantities at the temperature of interest T = T%. 
The average over the M a ^ partial values give the final results. In fact, we further improve 
the calculations and estimate the statistical deviations by performing this procedure (after 
relaxation) M rep times, so that in total the number of (a)-(b) MC steps is M tot = M eq + 
M afi x M rep . 

In (a), for each replica (at a distinct T n ), a site lattice / is chosen randomly. Then, its 
occupation variable o\ may change to a new value a[ according to the Metropolis prescription 
P = min{l, exp[— /3A"H]]^, where A% = H(cr') — H(cr) is the energy variation due to the 
occupation change. This is done until a full lattice covering and the process is repeated all 
over again M times, (b) In the second part, arbitrary pairs of replicas (say, at T n > and T n » 
and with microscopy configurations a' and a") can undergo temperatures switchings, with 
probability (f3 n = (A; B T n ) _1 ) 

p n W = min{l, exp[(/3 n , - M(n(a') - H(a"))}}. (2) 

The PT algorithm is schematic represented in Fig. 1 (a). 

Although the above prescription is rather simple, few technical aspects should be ob- 
served. First, it is necessary to find a good compromise between the p's values (which in- 
crease with AT /N decreasing) and the replicas number N. This is so to guarantee relatively 
frequent exchanges, while keeping the computational efforts low. Hence, extra procedures 
have been proposecr^i"— . Here we use only the ones explained above. However we men- 
tion that for our present systems, one of us has tested some of these extra implementations 4 
(always assuming arbitrary n"s and n"'s for the step (b) above), not finding any significant 
difference. Second, the system size (L) also imposes restrictions on the iV's. For small sys- 
tems, a few number of replicas is enough to assure rapid convergence. On the other hand, 
by increasing L the exchange probabilities (Eq. (j5j)) decreases, so the inclusion of extra 
copies becomes necessary. Such care has been explicit taken in our simulations. Finally, we 
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FIG. 1. Schematics of the (a) PT and (b) ST implementations. In this example, there have been 
two temperatures exchanges for the PT (Ti o T2 and T3 f> T5) and one temperature change for 
the ST (Ti -tTj). 



observe that most works that use the PT method implement the switching attempts only 
between adjacent replicas (i.e., at T n i and T n // =n / +1 ), in principle because the probability of 
exchanges decreases for increasing T n » — T n i. Nevertheless, it has been shown^ that non- 
adjacent exchanges are essential to speed up the crossing of high free-energy barriers (what 
we discuss in more details in Section II. C). Therefore, here we will allow exchanges between 
first (5 = 1), second (5 = 2), etc, neighbor replicas, meaning those between T n and T n+ $. 
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B. Simulated Tempering 



For the ST, a single realization of the model is considered, however, during the dynamics 
its temperature can assume the different values T n 's. The implementation is similar to that 
for the PT in Section II. A, but applied only to one copy of the system. Therefore, the 
previous step (b) now reads: A change T n > — > T n » may take place for the system according 
to the probability (with o its configuration) 

Pn'->n" = min{l, exp[(/3 n / - (3 n »)Ti(a) + (g n » - g n '))}- (3) 

The ST algorithm is illustrated in Fig. 1 (b). 

Note that p n '^ n " depends on the weights g's. Moreover, for a better sampling, the 
evolution should uniformly visit all the established temperatures. This is just the case when 
9n — Pn fm with f n the system free energy at T n — dlilk. To obtain / is not an easy task. For 
instance, in Ref.— its exact (numerical) values follows from /„ = — ln[Z n ]/(Vf3 n ), with the 
partition function Z n computed by an involving recursive procedure. Here, V is the system 
volume, which in a regular square lattice reads V = L 2 . In our examples we will consider 
this same protocol, but using a simpler numerical implimentation for Z n . Indeed, in the 
thermodynamic limit 

Z n = -(Ai°)) L , (4) 

where A^ ) is the largest eigenvalue of the transfer matrix T at T n (for details see, e.g., 
Ref.—). By its turn, A^ = (T(Sk, Sk)) / (5s k ,s k+1 ) can be calculated from straightforward 
Monte Carlo simulations^, where Sk is the lattice fc-layer configuration a\ t k, o%k, ■ ■ ■ , &L,k 
and Ss k ,s k+1 — 1 (= 0) if the k and k + 1 layers are equal (different). A central point is 
that in principle Eq. (pTJ) would hold true only for infinite size systems. However, if L is not 
too small, the above relation is extremely accurate and for any practical purpose gives the 
correct Z n , as we show in the next Section. Such way to determine p n '->n" will be named 
the ST (exact) free-energy method, ST-FEM. 

Finally, we observe that approximations for g are equally possible. One implementation 
being^ 

9n+l -9n~ (/W - P n )(U n+1 + U n )/2, (5) 

with U n = (Tin) ( n — 1; 2, . . . , N) the average energy at T n . The Ws can be evaluated from 
direct auxiliary simulations. For completeness we will also consider this ST approximated 
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(a) ST 




FIG. 2. Schematic illustration of the trajectories - succession of u's - generated by the PT and 
ST algorithms in the case of a complex topography for the relevant microstate space (resulting 
from specific parameters values). A higher sinuosity (usually associated to smaller T's) represents 
a higher difficulty to leave the particular region of S, full of energetic valleys and hills. The length 
of the paths are proportional to the number of steps taken by the algorithms. 

method, which we call ST-AM. 

C. The PT and ST methods near phase transition regimes 

The sampling of a statistical system when the phase space has a complicated landscape full 
of free-energy valleys and hills 30 is particularly delicate: one needs to uniformly visit different 
regions of (those more important for the given parameters), but which are separated 
by many entropic barriers^. In this case, the particular way in which a method evolves 
throughout the microstates space to generate Q - even with the use of enhanced procedures - 
may crucially determine the final outcome of sampling. For instance, non-ergodic "probing" 
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of the multiple domains^ can prevent the proper relaxation to equilibrium. 

The previous comments fit perfectly well first-order phase transitions, where the min- 
ima of the free-energy are separated by large barriers. Nevertheless, we observe that for 
second-order phase transitions, the divergence of time and spatial length correlations cre- 
ates strongly correlated configurations^. It leads to a certain clusterization of relevant parts 
of S at the critical point, with independent and unbiased Q difficult to obtain. So, although 
associated to different mechanisms, near both first and second order transitions we can ex- 
pect a "fragmented" phase space. Hence, even if the PT and ST are not crucially distinct 
in usual situations (in fact, the ST being slight better than the PT in few instances^), here 
we argue qualitatively that in such cases the PT can outperform the ST. 

Thus, for the above contexts of multiple basins^, the Fig. 2 schematically represents 
"stretches" of typical dynamical paths generated by the ST and PT algorithms. The suc- 
cessively visited oj's until leaving the domain - delimited by high local free-energy barriers 
(or cluster walls) - can form a very sinuous trajectory on that particular region of S due to 
a complex topography. 

Thus, consider first the ST, Fig. 2 (a). The initial microstate uq evolves (at T = 7\) in a 
very tortuous path, but in average towards the border of the domain, reaching u a after M 
steps. Then, it undergoes a temperature change 7\ — > T 3 and again evolves M steps getting 
to u' a , this time in a more straight trajectory because the higher T (note if there was no 
temperature change, the path would follow the dashed line displayed in the plot). Finally, 
there is a second successful attempt to change T, T 3 — > Tj > T 3 , and after M steps the 
system ends up very close to the barrier separating the basins. 

In Fig. 2 (b) we observe the PT dynamics, where just one successful temperature exchange 
takes place (between the only two replicas depicted). The microstate u>b iusd) is obtained 
from ojq after 2M steps at T = 7\ (T = Tj). Obviously, u' a in the ST must be in average 
closer to (farther from) the domain border than cj& (u^) in the PT implementation. Then, 
there is an exchange of temperatures and the evolution of Ud at T 1; after m < M steps, 
already makes the replica to cross the basin barrier to the microstate u' d . Furthermore, after 
At = M the state u b at Tj leads to a u' b close to the border. 

The above illustration - although certainly not extinguishing all the possibilities - is 
already representative of why the PT can be more efficient in sampling a space full of 
energetic valleys and hills (e.g., at phase transition regimes). It is so for the following 
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reasons: (i) In the PT, the existence of replicas at all the interval AT of temperatures 
generate some paths which more quickly will approach the domain borders, as seen in Fig. 
2 (b) for uq — > Ud at Tj. Moreover, the microstates along such trajectories at higher T's of 
course are usually more energetic, (ii) So, when finally there is an exchange of temperature, 
a microstate of high energy, even if now at lower T's, will demand a smaller number of steps 
to cross a barrier (like Ud — > oj' d in Fig. 2 (b)), and thus to start visiting other basins. On the 
other hand, trajectories of microstates of low energy, that during a certain At have evolved 
under small values of T's, e.g. ujq — > ujb in Fig. 2 (b), when shifting to higher temperatures 
will speed up their ways towards the barrier (ub — > u' h ). Note, nevertheless, that this is 
possible only if non-adjacent exchanges are allowed, the case we are assuming here, (iii) The 
above collective dynamics makes possible many of the replicas successfully leave a domain 
after fairly similar number of steps. Hence, once in another basin region, this "parallel" 
process can proceed in the same fashion, (iv) By its turn, we can face the ST as a "serial" 
process, then a faster drift towards the domain walls takes place only when T increases. 
As a consequence, the eventual more frequent temperature exchange for the ST— ~— not 
necessarily constitutes an advantage in complex S landscapes (as illustrated in Fig. 2). (v) 
Lastly, a not critical issue but which also may give some small advantage for the PT over 
the ST is that in the former, often the replicas (even at smaller T's) cross the domain high 
barriers more or less at the same time. Thus, once leaving a certain basin we already have 
a sample of microstates at T\ to make averages for the PT. As displayed in the Fig. 2 (a), 
for the ST it may happen that when the system reaches a microstate configuration able to 
cross the barrier, it is not at T\. Hence, an extra time is necessary for the system (naturally 
from the algorithm dynamics) to come back to Ti and so the averages to be performed. 

We finally observe that when the relevant space is more homogeneous in energy (e.g., 
far away from phase transitions), one should not expect so high increase of the trajectories 
sinuosity as we diminish T. Then, it is not difficult to realize that the listed differences 
between the PT and ST methods might not be important. 

The previous discussion is based on qualitative arguments. Of course, they should be 
corroborated by concrete quantitative studies. Next we analyze two systems near phase 
transition conditions. We will explicit show through detailed numerical simulations that 
indeed the PT algorithm is more efficient, specially in the case of first order phase transitions. 

10 



FIG. 3. For the Ising model with H = 0, L = 32, and units of J/ks, comparison between the 
partition function versus T calculated exactly^ and from Eq. (jl]). 
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FIG. 4. For the Ising model at T c , the auto-correlation functions versus r (in MC steps unities), 
simulated from the PT (continuous), ST-FEM (dashed) and ST- AM (dotted). 

III. THE ISING MODEL 



The model is defined by the following Hamiltonian 



v 



H 



(6) 



<i,j> i=l 

where < i, j > denotes nearest-neighbors pairs % and j of a d- dimensional lattice of V = L d 
sites. At each site i, the spin variable assumes the values cr, = ±1. J is the interaction 
energy and H is the magnetic field. The Ising model displays a second-order phase transition 
(ferromagnetic-paramagnetic) at T c ps 2.269 and H = 0. For a square lattice (d = 2), the 
transfer matrix diagonal elements are 

L 



T{S k , S k ) = exp 



(7) 



Our interest are in the energy u = (H)/V and modulus of the magnetization (which is 
the order parameter) m — (\ Yn=i a i\) /V P er volume. For their auto- correlation functions, 
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X (MC steps) x (MC steps) 

FIG. 5. For the Ising model at T c , the time evolution of u and m from a non- typical initial 
configuration simulated by the PT (continuous), ST-FEM (dashed) and ST- AM (dotted). 

we just set w = u and w = m in Eq. (JTJ. Regarding the parameters, we choose H = and a 
square lattice of L = 32. All the results are given in units of J /kg. To test the accuracy of 
the transfer matrix largest eigenvalue method in obtaing Z, in Fig. 3 we compare the exact 
partition function (obtained from the solution in Ref.— ) with that calculated from Eq. (J4]) for 
the Ising model and the above parameters. The agreement is indeed remarkable, indicating 
that even for L = 32, Z and consequently / is already very close to the thermodynamic 
limit value. 

Figure 4 displays C m and C u for T\ = T c . In the simulations we use only two replicas (with 
T2 = 2.4) and M — 1. From the plots we see that the auto-correlations decay faster when 
calculated by the PT than by both the ST- AM and ST-FEM methods. In Fig. 5 we compare 
the time evolution of the thermodynamic quantities starting from a "hard" initial condition, 
i.e., a configuration very different from the ones representative of the steady state. Thus, 
we consider a fully ordered configuration, which obviously is not typical at T = T c . This 
is a way of testing how efficient is a certain approach to drive the system to the stationary 
state. The Ising model at the transition temperature evolves to the equilibrium basically in 
the same fashion either when simulated by the PT or by both the ST's. 

So, we have that for a continuous phase transition (at least for the Ising model) the 
performances of the two tempering methods are essentially equivalent. Although at T c the 
PT shows faster auto-correlation decays (in contrast with the results of Ref for the same 
model, however calculated far away from the critical temperature), the stationary state is 
characterized by equivalent values of m and u for all methods. 
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IV. THE LATTICE-GAS MODEL WITH VACANCIES (BEG) 



A. Model 

The lattice-gas model (of size V = L d ) with vacancies is characterized by the Hamiltonian 

U = - E e r,s N r>i N s>j - E E ^ N r,i- (8) 

<i,j> r,s r i 

Here, r and s run over the species labels A and B, the e rs 's are the coupling energies (e AA , 
£bb, £ab and e BA ), N r>i = 0, 1 is the occupation numbers at site % for species r, and /i r 
is the species r chemical potential. The above model is equivalent to the Blume-Emery- 
Griffiths (BEG) spin-1 H— . Indeed, defining (with Oi = 0, ±1 the possible values for the 
spin variable) 

N Aii = {a} + a t )/2, N B , = {a} - a,)/2, (9) 

associating <Ji = 1 (-1) with the species A (B) and <7j = with a vacancy, and setting 
£aa = £bb and e AB = £ba, we get the BEG Hamiltonian 

« = - y E(J<Ti<T j + Ka*o?)- y E(Ha i -Do?), (10) 

<i,j> i 

for 

H = (ji A -fi B )/2, D = -(ji A + fi B )/2, 

J = (e AA - e AB )/2, K = {e AA + e AB )/2. (11) 

We will consider a square lattice with periodic boundary conditions. In this case, the 
transfer matrix diagonal elements read 

L 



T(S k , S k ) = exp [/? (( H + J vi+i,k) &i,k 



i=i 

+(J-D + K(l + af +1>k ))al k )]. (12) 

The model has two order parameters, q and m, defined by q = (J2j=i{^A,i + N B ,i))/V 
and m = (Y.Y=i(^A,i — N B j))/V. Also important is the quantity energy per volume, given 
by it = (H) i IV '. The auto-correlation are then obtained from w = q, w = m and w = u in 
Eq. (HJ. 
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B. Results 



For fixed K/J, H and T, the characteristic of the phase space is determined by D. In 
the regime we are interested, there are two phases if D is small, one rich in species A and 
the other in species B. For high values of D, the model displays a single gas phase, rich 
in vacancies. A strong first-order phase transition between these two situations takes place 
at D = D*, which obviously depends on K/J, H and T. For definiteness, in the following 
we study the BEG Hamiltonian assuming K/J = 3, H = and T = T\ = 1.4 (for other 
parameter values, see Sec. V). In such case, D* = 8.000 in the thermodynamic limit 4 . All 
the results will be presented in units of J/ks- 

It is well known that for different lattice-gas systems, approaches based on cluster 
algorithms^ are very appropriate to deal with metastability arising in first-order phase 
transitions. So, next we will compare results obtained from both tempering methods with 
those available from cluster calculations^. Regarding the parameters values, unless oth- 
erwise explicit mentioned, in the simulations we consider L = 20, D = 8.000 and the 
replicas in the temperature interval AT = 0.6. Also, whenever necessary we perform in 
total up to M tot = 8 x 10 7 simulation steps (see Sec. II. A) to evaluate the sought quantities. 
Furthermore, we always use M — 1. 

As the first comparative analysis, in Fig. 5 we plot the order parameter q probability 
distribution histogram for a long simulation run of 10 MC steps. As the chemical potential 
we set D = 8.004, instead of D = 8.000, since it leads to a same high for the two peaks of the 
bimodal order parameter probability distribution (we mention, nevertheless, that D = 8.000 
gives the same qualitative results). The agreement of the two tempering with the cluster 
method 6 is similar (in fact, a little better for the PT case). Such calculations show that 
for a long enough time, both the PT and ST are able to circumvent the metastable states, 
allowing the system to cross the free-energy barriers separating the different phases at the 
coexistence. 

Despite the previous agreement, the PT and ST do present differences when other aspects 
are analyzed. For instance, we show in Fig. 6 the time evolution of q towards the steady 
state, starting from a fully random initial configuration. We also consider distinct number 
of replicas N and temperature intervals AT. We find that under the same simulation con- 
ditions, generally the PT converges faster, being closer to the cluster results than the ST 
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FIG. 6. For the BEG model, the histograms of the order parameter q from a long simulation using 
the PT, ST and cluster algorithms. The insets are blow-ups of the (a) low and (b) high densities 
regions, q pa and q pa 1, respectively. 
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FIG. 7. For the BEG model, the time evolution of q from a fully random initial configuration, 
simulated from the PT, ST, and cluster. N denotes the number of replicas and AT = 0.6 if not 
otherwise specified in the curves. 

(ST-FEM and ST-AM). However, for the lower value of AT = 0.25, in all cases the system 
(up to 10 4 MC steps) cannot even escape the region near the initial random configuration. 
On the other hand, by increasing AT = 0.6 - although the probability for temperature ex- 
changing decreases - the system starts to move towards the stationary regime. Furthermore, 
the larger the number of replicas N, the faster the convergence. Finally we mention that 
the steady value of q = 2/3 at D = D* = 8.000 can be understood recalling that at the 
phase coexistence, two liquid phases (q ~ 1) coexist with one gas phase (q ~ 0). Since their 
weights are equal (1/3), we have q ~ 2/3 for any system size. 

Another interesting test is to perform the numerical simulations when the system is 
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FIG. 8. For the BEG model, m versus t in two distinct time intervals at the steady state (after 
M eq ), calculated with the PT, ST-FEM and ST-AM algorithm. 
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FIG. 9. Similar to Fig. 7, but comparing PT and cluster. 



already at the steady state. In Fig. 7 we show the time evolution of the "magnetization" 
m for both tempering methods at the phase coexistence. In the plots the time is shifted so 
to discard the M eq initial MC steps necessary for equilibration. We see that the tunneling 
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t x 10" 5 (MC steps) 

FIG. 10. For the BEG model, q versus the chemical potential, simulated from the PT (square), 
ST-FEM (triangle), ST- AM (circle), and cluster (x). The averages are taken at each 10 4 MC steps. 
In the inset, exactly the same curve but for the averages at each 5 x 10 4 MC steps. 

between the three different phases is substantially more frequent for the PT than for the ST. 
It being true along the whole evolution, as we have checked for an interval of 10 7 MC steps 
(in the Fig. 7 we show only two distinct simulation stretches). Actually, the PT tunneling 
pattern presents the same behavior than that observed in the notorious accurate cluster 
algorithm 5 , Fig. 8. Such results concrete exemplify some of the qualitative arguments given 
in Sec. II. C to explain why the PT should be more efficient than the ST around first-order 
phase transitions. 

A different efficiency for the methods is observed not just at the phase coexistence, but 
also for other values of the chemical potential D around D*. Figure 9 plots the order 
parameter q versus D for the PT and ST implementations, evaluating the averages at each 
M a ^ = 10 4 MC steps. Note that overall the PT is already quite close to the values obtained 
from the cluster algorithm, whereas both ST still show some discrepancy, specially for D > 
D*. If now the averages are calculate each M a j, = 5 x 10 4 MC steps, the ST also becomes 
closer to the cluster's (inset of Fig. 9). Once more such results can be understood in terms 
of the tunneling between the phases. For D ~ D*, we still can expect high free energy 
barriers. With the ST, the system does not cross such barriers a sufficient number of times 
if M a ^ = 10 4 . By increasing the number of MC steps for the averages, we generate a more 
representative Q and thus a better estimation for m. 

As a last efficiency measure, we consider the two relevant auto-correlation functions, 
C q (r) and C u (r), shown in Fig. 10. We should note that although time displaced correlation 
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FIG. 11. Auto-correlation functions versus r from the PT (continuous), ST-FEM (dashed) and 
ST- AM (dotted). 




FIG. 12. Mean probability of exchange versus the temperature T = T\ for the PT and ST-FEM. 
The symbols 5 = 1,2,3 refers, respectively, to exchanges allowed between first, second and third 
neighbors (see main text). 

functions are more commonly studied in the context of continuous phase transitions, in the 
present case they are an interesting auxiliary tool to compare the PT and ST performances. 
As it should be, the ST-FEM uncorrelates faster than the ST- AM. Nevertheless, we see that 
the C's decay even faster for the PT method (in fact, with a very drastic difference in the 
case of C u (t)). 

Usually, the frequency (measured in terms of a probability p*) in which a given tempering 
method changes the system temperature is taken as a good indication of its efficiency. For the 
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PT and ST algorithms, such quantity respectively reads 36 p* = (min{l, exp[(/?j — Pj)(?i(<Ti) — 
"H(ctj)]}) and p* = (min{l, exp[($ — /3j)"H(cr) + gj — gi}})- The averages are over T\, . . . , T/v, 
such that p* of order 5 is the mean from all the exchanges among T n and T n+ s (see Sec. 
II.A). 

In Fig. 11 we display p* as function of T = T x for the PT and ST-FEM (the ST- AM 
being similar to the latter), with A = 12 and AT = 0.55. As it can be seen, for any 5 the ST 
always presents a higher probability of acceptance than the PT, in agreement with previous 
studies^ 1 ^. Such findings are in contrast with our results here. Indeed, larger p*'s do not 
translate into a better performance of the ST, at least in the case of phase transitions as 
argued in Sec. II-C. Therefore, exchange probabilities alone should be faced with care when 
trying to characterize the best tempering method for a certain context. 

Finally, we show in Figs. 12 and 13 finite size analysis for the total density q and the 
isothermal susceptibility xt = /3L 2 ({q 2 ) — (q) 2 ) from the PT and ST-FEM. Continuous lines 
correspond to fitting curves by a method proposed in Ref.-. At the phase coexistence, ther- 
modynamic quantities scale with the system volume 3 ^* 3 ^. A discontinuous phase transition 
is characterized by a jump in the order parameter or even a delta function-like singularity for 
the susceptibility or specific heat. But this is so only at the thermodynamic limit. For finite 
systems not only the order parameter, but also other quantities are described by continuous 
functions^^. We should emphasizes that smooth curves are obtained only when one uses 
a simulation dynamics which correctly yields an appropriate sampling. For instance, from 
simple Metropolis algorithms, neither the crossing among isotherms nor accurate finite size 
analysis for smooth curves are possible. It is due to the presence of hysteresis effects^- 6 -, 
which hence demand tempering enhanced algorithm. From the plots we see that both the 
PT and ST give fairly good results. However, the cluster continuous curved is smoother and 
better fitted in the PT case, specially for the larger L = 30 value. 



V. REMARKS AND CONCLUSION 

In this paper we have presented a comparative study between two important enhanced 
sampling methods, namely, simulated (ST) and parallel (PT) tempering, considering spin- 
lattice models at phase transition conditions. Special attention has been payed to first-order 
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FIG. 13. q versus D for L equal to 10 (circle), 20 (square) and 30 (triangle), calculated from the 
(a) PT and (b) ST-FEM. Continuous lines are fitting results^. The curves collapse if plotted as 
qx (D- D*)L 2 (insets). 




FIG. 14. Susceptibility versus D for L equal to 10 (circle), 20 (square) and 30 (triangle), calculated 
from the (a) PT and (b) ST-FEM. The curves collapse if plotted as xt/L 2 x (D - D*)L 2 (insets). 



phase transitions at low temperatures (for the BEG model). In such regimes, more standard 
algorithms often give poor results because their difficulties to overcome the large free-energy 
barriers in the phase space, leading, e.g., to ergodicity breaking and artificial algorithm- 
induced hysteresis. We also have investigated the less critical case of second order-phase 
transition - for which no free-energy barriers exist but there is the formation of strongly 
correlated clusters (basin regions) 33 - for the well understood Ising model. 

As for the tempering implementations, we have followed the usual PT procedure, but 
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allowing temperature exchanges between non-adjacent replicas. For the temperature change 
probability weights g in the ST, we have assumed a recent proposed approximation 1 : 2 (ST- 
AM) and a new alternative exact approach (ST-FEM), based on the eigenvalues of the 
transfer matrix 2 ^. The ST-FEM here is formally similar to that in Ref.— , but avoids the 
necessity to implement more complicated recursive procedures to estimate the partition 
function. 

Different comparative analysis, both at the transient regime and already at the steady 
state, have been carried out. Despite the facts that: (i) after long times (thus demanding 
large computational effort) the final results from the PT and ST are similar; and (ii) the PT 
displays a smaller exchange probability than the ST; we have found that for discontinuous 
phase transitions the PT is always more efficient in any verified aspect. The main reason for 
this is basically that the PT enables the system to cross free-energy barriers more frequently 
than the ST: either at or near phase coexistence conditions (as explicit illustrated, e.g., in 
Figs. (7) and (8)). Furthermore, besides the quantitative numerical results, we also have 
presented heuristic arguments for why it should be expected. 

Results for the instructive Ising model at the critical temperature (second-order phase 
transition) have also agreed with our qualitative predictions. Indeed, far away from T c it 
has been reported a faster convergence for the ST—. We have shown that for T ~ T c just 
the opposite takes place, with the auto-correlations decaying faster for the PT. 

For completeness, we also have analyzed other values of Kj J for the BEG model (not 
shown), in particular for Kj J = 0, the so called Blume-Capel model. The calculations at 
the first-order transition (7\ = 0.4 and D = 1.9968) have corroborated the higher efficiency 
of the PT over the ST. More specifically, until M tot = 3 x 10 7 , the system when simulated 
with the ST-AM has not reached the steady state, whose values for the thermodynamic 
quantities were different from those obtained by the ST-FEM, PT and cluster algorithms. 
Furthermore, the ST-FEM have agreed with the PT and cluster only for long M to tS. Time- 
displaced correlation functions decays and actual thermodynamic quantities convergence 
were always faster for the PT. 

A second contribution of this work has been an (numerically simpler) alternative way to 
calculate the exact g in the ST method. When comparing the ST-AM with the ST-FEM, 
we have found that the ST-FEM allows the system to converge to steady regime quicker 
than the ST-AM (see above). In addition, at the steady state, configurations generated by 
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ST-FEM uncorrelate faster than those by the ST-AM. On the other hand, with respect to 
the frequency in which the system tunnels between different phases at the coexistence and 
the final sough thermodynamic quantities, both implementations are similar, but the latter 
only for long M tot 's. 

Summarizing, at phase transition regimes the PT and ST provide the same results for 
long (sometimes even costly) simulations. However, we find that for all the tested measures, 
the parallel converges faster than the simulated tempering. Also, even in such situation of a 
better performance from the PT, still the rate of temperature switching is higher for the ST. 
Thus, another message from our work is that alone, the switching rates are not sufficient to 
characterize the efficiency of a tempering enhanced sampling algorithm. 
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